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This  research  program  investigated  signal  processing  problems  encountered  in  high- 
resolution  image  formation.  Reliable  imaging  of  scenes  with  high  resolution  and  high  speed 
is  an  important  and  key  part  of  any  defense  system.  In  some  imaging  systems,  the  image 
has  to  be  constructed  from  linear  measurements  and  convex  constraints  (such  as  upper  and 
lower  bounds  on  image  sample  magnitude,  and  support  limits).  We  studied  iterative,  finite 
parameter  reconstructions  that  lead  to  images  that  meet  the  constraints,  match  the  data  to 
within  a  pre-specified  tolerance,  and  come  closest  to  a  given  nominal. 

We  addressed  three  different  aspects  of  signal/image  reconstruction,  namely: 

1.  Developing  fast  algorithms  for  high  resolution  signal/image  reconstruction. 

2.  Resolutions  Analysis  of  signal  reconstruction  algorithms. 

3.  Numerical  and  computational  aspects  of  signal  reconstruction  viz.  stability  and  regu¬ 
larization. 
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Important  results  in  each  of  the  above  three  topics  are  summarized  below. 

Summary  of  Important  Results 

1.  Fast  Algorithms 

The  standard  approach  to  the  problem  of  image  reconstruction  from  linear  measurements 
and  convex  constraints  has  been  an  alternating  projections  onto  convex  sets  (POCS)  algo¬ 
rithm,  studied  by  Youla  and  others.  The  approach  suffers  from  slow  (linear)  convergence, 
high  computational  cost  and  non-unique  limits.  We  have  developed  a  quadratically  con¬ 
vergent  iterative  algorithm  (Newton  algorithm)  for  this  problem.  Central  to  the  Newton 
algorithm  is  the  derivative  of  the  nonlinear  projection  operator  onto  a  convex  set.  We  ob¬ 
tained  a  new  general  mathematical  result  for  the  existence  and  construction  of  the  derivative 
of  the  projection  operator  for  a  class  of  convex  sets.  This  result  was  then  used  to  give  the 
Newton  algorithm  for  the  signal  recovery  problem.  A  salient  feature  of  the  .algorithm  is  the 
quadratic  rate  of  convergence. 

Avallabllttir 

UySiT 
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We  also  studied  the  implementation  aspects  of  the  algorithm  and  developed  a  computation 
and  memory  efficient  implementation  of  the  algorithm  using  conjugate-gradient  iterations 
within  each  Newton  iteration.  A  remarkable  feature  of  the  algorithm  with  this  implemen¬ 
tation,  is  that  each  iteration  has  similar  computational  complexity  as  an  iteration  of  the 
POCS  algorithm.  The  faster  rate  of  convergence  of  the  algorithm  (as  compared  to  POCS) 
thus  enables  us  to  compute  high  resolution  reconstructions  with  fewer  computations.  The 
algorithm  has  been  tested  extensively  on  imaging  examples  and  we  have  obtained  extremely 
good  performance.  In  many  image  formation  schemes,  6  to  7  iterations  seem  to  suffice, 
compared  to  the  100s  of  iterations  required  by  classical  methods  such  as  POCS. 

2.  Resolution  analysis 

Resolution  ability  is  the  ability  to  reproduce  fine  details  such  as,  narrow  peaks  or  closely 
spaced  peaks  in  a  signal.  The  study  of  resolution  is  important,  since  understanding  the 
relationship  between  resolution  limits  and  the  various  components  of  a  recovery  problem 
and  algorithm,  could  help  us  design  better  data  acquisition  schemes  and  algorithms. 

The  earliest  definition  of  resolution  limit  is  the  Rayleigh  Resolution  Limit.  This  definition 
is  based  solely  on  the  observed  data  and  not  on  any  recovery  algorithm.  The  definition  is 
acceptable  when  there  is  no  processing  of  the  data  to  recover  or  enhance  the  features  based 
on  exploiting  prior  information.  The  Rayleigh  limit  is  thus  a  lower  bound  on  the  achievable 
resolution. 

We  find  that  where  infinitely  many  noise-free  measurements  are  available,  the  resolution 
achievable  is  in  fact  independent  of  the  width  of  the  sampling  pulse  and  depends  only  the 
inter-sample  distance.  The  Rayleigh  limit,  on  the  other  hand,  is  dictated  by  the  width  of  the 
sampling  pulse.  In  the  presence  of  observation  noise,  however,  the  notion  of  exact  recovery 
has  to  be  abandoned  and  a  new  measure  of  resolution  is  necessary.  We  define  a  new  measure 
based  on  allowable  levels  of  worst-case  error,  and  find  that  the  resolution  limit  depends  on 
the  method  of  regularization  used  in  the  recovery  algorithm. 

In  the  more  practical  situation  in  which  only  finitely  many  noisy  observations  are  avail- 
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able,  the  worst-case  error  is  unbounded  and  so  we  have  to  restrict  the  search  to  a  smaller 
set  of  signals.  In  order  to  meaningfully  describe  resolution  limits  for  the  problem  of  signal 
recovery  from  finitely  many  noisy  observations,  we  restrict  the  class  of  signals  to  the  set  of 
bandlimited  and  essentially  timelimited  signals,  since  it  describes  most  signals  encountered  in 
practice.  This  set  is  characterized  by  the  well  known  orthonormal  family  of  functions  called 
the  Prolate  Spheroidal  Wave  Functions  and  is  known  to  be  approximately  finite  dimensional, 
which  enables  us  to  seek  reconstructions  from  a  lower  dimensional  subspace  of  the  space  of 
bandlimited  signals.  Reduction  in  dimension  causes  an  error  in  the  reconstruction,  which  we 
call  the  intrinsic  error.  A  second  error  is  incurred  while  determining  the  parameters  describ¬ 
ing  the  lower  dimensional  reconstruction.  The  reconstruction  error  is  then  the  sum  of  these 
two  errors.  We  show  that  the  worst-case  values  of  these  two  errors  can  be  pre-computed  for 
each  choice  of  reduced  dimension.  The  error  computation  provides  both  an  optimal  choice 
of  dimension  and  a  precomputed  bound  on  the  resolution  ability  of  the  algorithm. 

3.  Numerical  and  Computational  aspects: 

Most  reconstruction  problems  are  ill-posed.  Therefore,  regularization  is  needed  for  an 
approximate  but  stable  computation  of  the  solution.  The  truncated  S  VD  approach  is  popular 
method  of  regularization.  The  selection  of  the  proper  rank  is  the  main  problem  in  this 
technique.  We  have  developed  a  new  rank  selection  strategy  based  on  a  bound  on  the  noise 
energy  alone.  In  the  absence  of  bounds  on  the  signal  energy  we  enforce  our  belief  that  the 
signal  energy  cannot  be  large,  by  introducing  a  penalty  on  the  solution  in  the  worst-case 
(min-max)  analysis.  The  rank  is  selected  so  that  the  worst  case  error  with  a  penalty  on 
the  solution  norm  is  minimum.  This  method  of  regularization  was  tested  in  bandlimited 
extrapolation  problems  and  it  gave  excellent  results. 
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Resolution  Limits  in  Signal  Recovery* 

S.  Dhaxamipragada^  and  K.  S.  Arun^ 
September  27,  1994 


Abstract 

Resolution  analysis  for  the  problem  of  signal  recovery  from  finitely  many  lineatr 
samples  is  the  subject  of  this  paper.  The  classical  Rayleigh  limit  serves  only  as  a  lower 
bound  on  resolution  since  it  does  not  assume  any  recovery  strategy  and  is  based  only  on 
observed  data.  We  show  that  details  finer  than  the  Rayleigh  limit  can  be  recovered  by 
simple  linear  processing  that  incorporates  prior  information.  We  first  define  a  measure 
of  resolution  based  on  allowable  levels  of  error  that  is  more  appropriate  for  current  signal 
recovery  strategies  than  the  Rayleigh  definition.  In  the  practical  situation  in  which  only 
finitely  many  noisy  observations  are  available,  we  have  to  restrict  the  class  of  signals 
in  order  to  make  the  resolution  measure' meaningful.  We  consider  the  set  of  bandlim- 
ited  and  essentially  timelimited  sign2Js  since  it  describes  most  signals  encountered  in 
practice.  For  this  set  we  show  how  to  precompute  resolution  limits  from  knowledge  of 
measurement  function2ils,  signal-to-noise  ratio,  passbemd,  energy  concentration  regions, 
energy  concentration  factor,  and  a  prescribed  level  of  error  tolerance.  In  the  process  we 
alro  derive  an  algorithm  for  high  resolution  signal  recovery.  We  illustrate  the  results 
with  examples  in  one  and  two  dimensions. 
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1  Introduction 

The  problem  of  recovering  signals  from  linear  measurements  arises  in  many  applications,  and 
several  algorithms,  linear  and  nonlinear,  have  been  developed  -and  analyzed  for  this  problem 
[1, 2,  3, 4, 5, 6, 7,  8,  9].  However,  the  fundamental  question  regarding  the  resolution  ability  of 
a  recovery  algorithm  is  often  left  unanswered.  Resolution  ability  is  the  ability  to  reproduce 
fine  details  such  as,  narrow  peaks  or  closely  spaced  peaks  in  a  signal.  The  study  of  resolution 
is  important;  in  many  applications  it  is  necessary  for  the  reconstruction  algorithm  to  have 
a  certain  minimum  resolution  in  order  to  be  elfective.  For  example,  consider  the  digital 
mammography  application,  where  a  2-D  X-ray  profile  is  to  be  reconstructed  from  samples. 
With  current  sensor  technology  and  physical  limitations,  the  sampling  operation  amounts  to 
2D  pulse  sampling  with  pulse  width  (A)  at  least  50  microns  and  sample  spacing  (r)  at  least 
50  microns,  which  limits  the  “resolution”  obtainable  from  a  single  exposure  to  50  microns 
or  larger,  i.e.,  details  in  the  image  that  are  narrower  than  50  microns  cannot  be  reproduced. 
The  integration  with  a  pulse  of  width  50  microns  causes  smearing  of  finer  details,  and  the 
sample  spacing  of  50  microns  can  cause  us  to  miss  these  details  altogether.  However,  early 
detection  of  breast  carcinoma  requires  that  features  of  width  25  microns  be  reproduced  in 
the  image.  The  grain  size  of  an  X-ray  film  is  small  enough  for  these  features  to  appear  in 
the  more  conventional  analog  X-ray  photographs.  To  make  digital  mammography  equally 
useful  for  diagnostic  radiology,  resolution  of  at  least  25  microns  is  required.  More  generally, 
the  study  of  resolution  limits  is  important  since,  it  could  help  us  assess  the  effectiveness  of 
a  particular  algorithm,  and  compare  different  algorithms  in  a  rational  manner.  Moreover, 
understanding  the  relationship  between  resolution  limits  and  the  various  components  of  a 
recovery  problem  and  algorithm,  could  help  us  design  better  data  acquisition  schemes  and 
algorithms. 

The  problem  of  resolution  analysis  is  twofold:  first,  we  need  a  meaningful  measure  of 
resolution  ability,  and  second,  we  have  to  be  able  to  analyze  the  performance  of  a  recon¬ 
struction  algorithm  in  terms  of  the  defined  resolution  measure.  The  earliest  definition  of 
resolution  limit  is  the  Rayleigh  Resolution  Limit.  It  is  defined  as  follows: 

Deffnition  1.1  (Rayleigh  Resolution  Limit)  [10]  If  two  equally  strong  point  sources 
(impulse  intensities) ,  5  or  more  apart,  are  reproduced  as  peaks  with  at  least  a  19%  intensity 
dip,  and  if  sources  less  than  5  apart  are  not  reproduced  as  well,  the  resolution  limit  is  said 
to  be  5.  • 

Equivalently,  the  Rayleigh  resolution  limit  is  inversely  proportional  to  the  main  lobe  width 
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of  the  point  spread  function  of  the  blurring  (or  sampling)  kernel.  This  definition  is  based 
solely  on  the  observed  data  and  not  on  any  recovery  algorithm.  The  definition  is  acceptable 
when  there  is  no  processing  of  the  data  to  recover  or  enhance  the  features  based  on  exploiting 
prior  information.  The  Rayleigh  limit  is  thus  a  lower  bound  on  the  achievable  resolution. 
We  might  be  able  to  do  better  with  clever  signal  processing  that  exploits  prior  information, 
but  we  should  always  be  able  to  achieve  at  least  as  much  resolution  as  specified  by  the 
Rayleigh  limit. 

We  find  that  where  infinitely  many  noise-free  measurements  are  available,  the  resolution 
achievable  is  in  fact  independent  of  the  width  of  the  sampling  pulse  and  depends  only  the 
inter-sample  distance.  The  Rayleigh  limit,  on  the  other  hand,  is  dictated  by  the  width 
of  the  sampling  pulse.  In  the  presence  of  observation  noise,  however,  the  notion  of  exact 
recovery  has  to  be  abandoned  and  a  new  measure  of  resolution  is  necessary.  We  define  a 
new  measure  based  on  allowable  levels  of  worst-case  error,  and  find  that  the  resolution  limit 
depends  on  the  method  of  regularization  used  in  the  recovery  algorithm. 

In  the  more  practical  situation  in  which  only  finitely  many  noisy  observations  are  avail¬ 
able,  the  worst-case  error  is  unbounded  and  so  we  have  to  restrict  the  search  to  a  smaller 
set  of  signals.  In  studying  the  problem  of  resolution  in  signal  recovery  Root  et  al.  [11],  [12] 
recognized  the  need  for  finite-dimensional  approximations  to  overcome  instability.  They  in¬ 
troduce  the  concept  of  an  error  number,  which  is  defined  as  the  mean-squared  error  averaged 
by  the  dimension  of  the  approximating  subspace.  As  dimension  increases,  the  “detail”  in 
the  estimate  increases  but  so  does  the  error  number.  The  authors  thus  bring  out  the  trade¬ 
off  between  the  achievable  level  of  detail  (i.e.,  resolution)  and  the  error  in  reconstruction. 
However  they  do  not  quantify  resolution  (detail)  and  fail  to  provide  resolution  bounds. 

We  appeal  to  the  Fourier  uncertainty  principle  to  bring  out  the  relationship  between 
resolution  (detail)  and  bandwidth.  In  this  sense  our  is  simil&r  in  spirit  to  the  classical 
Rayleigh  resolution  limit,  but  is  based  on  a  prescribed  tolerance  of  the  relative  error. 

In  order  to  meaningfully  describe  resolution  limits  for  the  problem  of  signal  recovery 
from  finitely  many  noisy  observations,  we  restrict  the  class  of  signals  to  the  set  of  bandlimited 
and  essentially  timelimited  signals,  since  it  describes  most  signals  encountered  in  practice. 
This  set  is  characterized  by  the  well  known  orthonormal  family  of  functions  called  the 
Prolate  Spheroidal  Wave  Functions  and  is  known  to  be  approximately  finite  dimensional, 
which  enables  us  to  seek  reconstructions  from  a  lower  dimensional  subspace  of  the  space  of 
bandlimited  signals.  Reduction  in  dimension  causes  an  error  in  the  reconstruction,  which 
we  call  the  intrinsic  error,  A  second  error  is  incurred  while  determining  the  parameters 
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describing  the  lower  dimensional  reconstruction.  The  reconstruction  error  is  then  the  sum 
of  these  two  errors.  We  show  that  the  worst-case  values  of  these  two  errors  can  be  pre¬ 
computed  for  each  choice  of  reduced  dimension.  The  error  computation  provides  both  an 
optimal  choice  of  dimension  and  a  precomputed  bound  on  the  resolution  ability  of  the 
algorithm. 

This  paper  is  organized  as  follows:  In  Section  2  we  formulate  the  signal  recovery  problem 
in  a  vector  space  setting.  In  Section  3  we  analyze  the  resolution  limits  for  the  ideal  situation 
of  infinitely  many  noise-free  observations  and  describe  the  effects  of  noise  and  regulariza¬ 
tion.  Based  on  the  results  obtained,  we  suggest  a  method  for  improving  resolution  in  the 
mammography  application.  In  Section  4,  which  is  the  bulk  of  the  paper,  we  examine  the 
practical  situation  of  finitely  many  noisy  observations  for  resolution  limits.  We  demonstrate 
the  results  in  one  and  two  dimensional  examples. 

2  The  Signal  Recovery  Problem 

We  consider  the  problem  of  reconstructing  ID  continuous-index  signals  from  discrete  linear 
measurements.  The  results  presented  here  can  be  easily  generalized  to  multidimensional 
signals.  We  demonstrate  the  generalization  by  way  of  an  example  in  Section  4. 

Let  L2{R)  be  the  space  of  finite-energy  continuous-index  signals  with  the  natural  inner 
product  defined  by 

JR 

where  the  overbars  denote  complex  conjugation.  Let  Bs  be  the  subspace  of  all  signals 
ha.nd1inriit.ftd  to  Ps  =  [^,f]  and  let  B  denote  the  orthogonal  projection  operator  onto  Bs- 
In  practical  terms,  B  is  simply  a  ideal  bandpass  filter  with  a  passband  Ps  =  [^,  f].  Let 
A’(n)  denote  the  continuous-time  Fourier  transform  (CTFT)  of  x{t), 

X(fi)  =  T{x{t)}  =  r  x{t)e-^^*dt, 

•/— OO 

and 

We  address  the  problem  of  recovering  a  signal  from  Bs  by  using  discrete  linear  measure¬ 
ments.  Every  linear  continuous  measurement  functional  on  L2{R)  can  be  expressed  as  an 
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inner  product  with  a  measurement  signal  in  L2{R)  [13].  Let  gi  be  measurement  signals 
giving  measurements  yd{i)  ^ 

yd{i)  =  gi)  =  L  x{t)gi{t)dt.  (1) 

J  R 

Let  T  be  the  linear  bounded  operator  on  Bs  representing  the  measurement  process. 
Then 

Tx  =  y,,  xe  Bs,  (2) 

where  yd  is  the  vector  of  measurements.  If  the  number  of  measurements,  p,  is  finite,  yd  lies 
in  CP,  otherwise  yd  lies  in  hiZ),  the  space  of  finite-energy  discrete-index  signals  with  the 
inner  product  _ 

{xd,  yd)h  =  13  ^d{n)yd{n).  (3) 

n&Z 

We  will  use  the  subscript  ‘d’  to  identify  discrete-index  signals,  their  transforms,  and  discrete¬ 
time  frequency  responses.  Let  Xd  €  L2([— Jr,  tt])  be  the  discrete  time  Fourier  transform  of 
Xd  €  hiZ), 

Xd{(^)  =  13  Xdin)e~^'^,  w  6  [-tt,  tt]. 

Finally,  the  adjoint  operator  T*  maps  a  vector  Vd  6  C^  or  l2{Z)  to  a  signal  T*Vd  = 

13  Vd{i)Bgi  in  Bs,  a  simple  linear  combination  of  the  measurement  signals.  Here  M  denotes 
M 

either  the  index  set  {1, 2,  •  •  • ,  p}  or  Z. 

In  the  ideal  situation  of  accurate  measurements,  the  problem  of  signal  recovery  from 
linear  measurements  is  equivalent  to  that  of  finding  a  solution  to  the  linear  operator  (2). 
However,  in  practice,  the  measurements  are  corrupted  by  noise.  Let  nj  denote  the  noise 
vector,  then 

Zd  =  Tx  +  rid  (4) 

and  the  signal  recovery  problem  is  that  of  reconstructing  x  £  Bs  from  Zd. 

In  the  next  section  we  consider  the  signal  recovery  problem  where  infinitely  many  mea¬ 
surements  are  available  and  examine  the  resolution  limit  of  the  minimum  norm  least  squares 
solution  to  (2).  We  also  describe  the  elfects  of  noise  on  the  resolution  limit  in  this  situation. 
The  bulk  of  this  paper,  however,  deals  with  resolution  limits  in  the  practical  case  of  a  finite 
number  of  noisy  measurements,  starting  from  Section  4 

3  Resolution  Limit  with  Infinitely  Many  Measurements 

It  is  well  known  that  if  all  signals  bandlimited  to  [“ji'f]  can  be  reconstructed  perfectly, 
then  two  point  sources  spaced  5  apart  will  show  up  as  distinct  peaks  in  the  reconstruction. 
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This  observation  is  a  natural  basis  for  our  first  definition  of  resolution. 

Definition  3.1  (Resolution  limit  under  ideal  conditions)  A  reconstruction  algorithm 
is  said  to  have  an  ideal  resolution  of  5  if  signals  bandlimited  to  [— j,  f  ]  can  be  reconstructed 
perfectly  under  noise-free  conditions.  * 

This  definition,  unlike  the  Rayleigh  definition,  is  based  on  the  recovered  signal  instead  of 

on  the  observed  signal.  Hence,  the  resolution  limit  will  depend  on  the  recovery  strategy 
adopted.  We  will  show  next  that  we  can  do  better  than  the  Rayleigh  resolution  limit. 

When  (2)  admits  no  exact  solution  because  of  noise  in  measurements,  a  popular  recourse 
is  to  seek  the  unique  minimum  norm  least  squares  solution  of  (2) 

£mnls  =  T*{TT*)^z,,  (5) 

where  (Tr*)t  denotes  the  pseudoinverse  of  the  composition  TT*.  From  the  description  of 
the  measurement  operator  and  its  adjoint  it  follows  that,  for  uj  € 

iTT-v,){i)  =  Yl  z4k){Bgk, gi)L,(Ry  (6) 

k&2 

In  many  applications  such  as  digital  mammography  and  deconvolution,  the  measurement 
functions  are  uniformly  translated  versions  of  a  basic  sampling  kernel/pulse  go: 

gk{t)  =  9oit  -  *r),  r  >  0.  (7) 

With  this  assumption,  {Bgk,gi)L2(R)  ~  {Bgk^  Bgi)i^(^R)  depends  only  on  i  —  k  and  not  on 
absolute  time,  and,  hence,  TT*  is  a  convolution.  The  shift-invariant  property  of  the  sam¬ 
pling  functions  and  the  TT*  operator  is  exploited  to  obtain  frequency  domain  expressions 
for  r,  T*,  TT*  and  the  MNLS  solution  next. 

From  (1)  and  (7)  we  obtain, 

(8) 

keZ 

where  Yi(w)  is  the  DTFT  of  yd  and  A’(fi)  and  <j(n)  are  the  CTFTs  of  x  and  go  respectively. 
If  x(t)  =  T*Vd,  then 

x(t)  =  {T*Vd)  {t)  =  x;  Vdik)  (Bgo)  it  -  kr)  (9) 

which  in  the  frequency  domain  is 

A(fi)  =  n(n)G(fi)F<i(Dr).  (10) 
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Here  n(fi)  is  the  frequency  response  of  the  ideal  bandpass  filter  for  passband  P5  =  ,  f], 

XCfi)  is  the  CTFT  of  x  and  Vi{u)  is  the  DTFT  of  Vd-  The  TT*  operation  is  equivalent  to 
discrete-time  convolution  with  a  kernel  given  by 

h{k)  —  {Bgk,Bgo)i^(^R)  =  f  go{t  —  kT){Bgo){t)dt.  (11) 

J  R 


Thus,  if  U(i  =  {TT*)zd  then 
where 


Vdiu)  =  Hd{uj)Zd{uj), 


keZ 


(12) 

(13) 


The  (TT*)^  operation  is  thus  equivalent  to  discrete- time  filtering  with  a  filter  of  frequency 
response 

f  1  ir./,  -X  n 

(14) 


{0^ 


„5 

Hi(u)  =  0 


Combining  the  frequency  domain  expressions  for  T*  and  (rr*)^,  we  find  that  the  minimum 

norm  least  squares  (MNLS)  solution  x  =  T*{TT*)^Zd  can  be  computed  in  the  frequency 

domain  as  ^  ^ 

X(fi)  =  Rin)Zd{nT),  (15) 


where 


where  ff(i(flr)  ^  0 
where  Hd{^T)  =  0 


It  is  remarkable  that  the  MNLS  algorithm  is  a  linear,  time-invariant  filter  even  in  the 
absence  of  any  restrictions/bounds  on  the  sampling  rate  (^)  or  pulse  width  A. 

Conventional  wisdom  holds  that  the  resolution  5  is  limited  by  both  r  and  A,  i.e.,  that  to 
achieve  a  resolution  of  5,  not  only  must  the  sample  spacing  t  be  less  than  S,  but  the  pulse 
width  A  must  also  be  smaller  than  S.  The  first  condition  is  considered  necessary,  because 
the  Nyquist  sampling  frequency  for  a  passband  [^,  j]  is  ^  (which  means  the  minimum 
value  of  sample  spacing  r  is  5).  The  second  condition  is  deemed  necessary  because  of  the 
Rayleigh  resolution  limitations.  Two  point  sources  (Dirac  delta  functions)  spaced 
S  apart  are  reproduced  in  the  measured  signal  yd  zs  a,  pair  of  superimposed  copies  of  the 
measurement  kernel  go  with  a  relative  translation  of  5.  If  the  width  A  of  kernel  go  is  less 
than  5,  the  two  copies  do  not  overlap  and  the  two  point  sources  stand  out  as  resolved 
distinct  peaks  in  yd- 

The  following  rather  simple  analysis  of  the  MNLS  reconstruction  (15)  shows  that  (at 
least  in  this  ideal  noiseless,  infinite-data  situation),  the  Rayleigh  limit  is  irrelevant,  and  that 
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even  if  the  two  point  sources  do  not  appear  to  be  resolved  in  the  measured  signal  yd,  they 
may  be  resolved  in  the  reconstruction  x.  The  analysis  shows  that  in  this  ideal  situation, 
resolution  is  determined  solely  by  the  sample  spacing  (i.e.,  5  =  t)  and  is  independent  of 
pulse  width  A  and  even  the  shape  of  the  sampling  kernel  qq. 

If  r  is  smaller  than  the  Nyquist  sampling  interval  for  Bs  =  ['X’f]  ^ 

dlimited  to  ^]),  then  there  will  be  no  aliasing  in  Hdit^)  in  (13).  In  this  case, 

Hi(u)  =  n(^)|G(^)|^ 

(16) 


Therefore, 


(  n(mG((i)  _  1 

R{Q)  =  I  n(n)tG(n)p  -  ^ 


when  n(f2)G(fi)  ^  0 
when  n(n)G(n)  =  0 


(17) 


and,  assuming  noise-free  observations, 

if(n)  =  { 

From  this  analysis  it  is  clear  that,  under  noise-free  conditions  and  infinitely  many  sam¬ 
ples,  the  reconstructed  spectrum  X(fi)  differs  from  the  true  spectrum,  X(fi)  only  at  the 
frequencies  at  which  G(fi)  =  0.  If  g{t)  is  of  finite  duration,  (A  <  oo),  the  Fourier  Uncer¬ 
tainty  Principle  states  that  G(fi)  cannot  be  bandlimited  i.e.  the  set  of  frequencies  where 
G{Q)  =  0,  has  zero  measure.  Thus  X(fi)  =  X(fi)  almost  everywhere  in  this  case,  which 
means  that  the  minimum  norm  reconstruction  algorithm  has  an  ideal  resolution  equal  to  r 
and  not  pulse  width  A.  This  result  is  independent  of  the  shape  of  the  sampling  pulse,  as 
long  as  its  support  A  is  finite. 

From  the  above  analysis,  we  now  see  that  it  is  possible  to  improve  resolution  in  digital 
mammography  and  other  applications  in  which  the  width  A  of  the  sampling  pulse  go  is 
greater  than  the  desired  resolution,  5,  by  the  use  of  multiple  exposures  based  on  staggering 
the  sampling  grid  by  r,  r  <  5,  with  every  new  exposure.  This  strategy  is  illustrated 
in  Figure  1  for  a  ID  problem  with  a  square  sampling  pulse.  For  simplicity  of  analysis, 
we  assume  that  the  sample-spacing  within  an  exposure  is  mr,  where  m  is  the  number  of 
exposures.  Thus,  after  the  m  exposures  are  interleaved,  we  obtain  uniform  sampling  with 


when  n(D)G(D)  ^  0 
when  n(f2)G(fi)  =  0 


(18) 
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Figure  1:  The  multi-exposure  sampling  strategy. 


spacing  t  which  can  be  a  small  fraction  of  A.  In  a  2D  application  like  digital  mammography, 
this  strategy  requires  multiple  exposures  with  2D  translations  of  the  collimators  (relative 
to  the  object)  between  exposures.  For  instance,  to  improve  resolution  from  50  microns  to 
25  microns  in  each  dimension  in  this  2D  problem,  we  will  need  4  staggered  exposures.  A 
similar  sampling  strategy  was  suggested  for  tomography  and  deconvolution  applications  in 
[14,  15]. 


3.1  Effect  of  noise  on  resolution 


In  practice  measurements  are  usually  corrupted  by  noise  and,  generally,  only  a  finite  number 
of  messurements  are  available,  which  causes  the  resolution  limit  to  be  affected  by  the  shape 
and  width  of  the  sampling  pulse  and  by  the  signal-to-noise  ratio  (SNR).  The  effects  of 
observation  noise  on  the  reconstruction  will  be  considered  next  while  we  still  assume  that 
infinitely  many  measurements  are  available. 

If  Nd{(^)  is  the  DTFT  of  n^,  the  noise  vector,  then  in  the  adequately  sampled  case  (i.e., 
T  <  5)  the  MNLS  solution  in  the  frequency  domain  is  given  by 


UsLsm  =  I 


where  n(n)G(f2)  ^  0 
where  n(D)G(f2)  =  0 


(19) 


When  the  magnitude  of  G(fi)  is  small,  the  noise  component  (n(nf|^n)p')  becomes  rela¬ 
tively  large,  which  leads  to  a  totally  unacceptable  and  unstable  solution.  Therefore,  regular¬ 
ization  is  required  for  an  approximate  but  stable  solution.  A  simple  regularization  scheme 
is  to  zero  the  solution  whenever  |G(f2)|  falls  below  a  certain  value.  The  regularized  MNLS 
solution  in  this  case  is 


=  I  n(SlG(n)|- 


where  in(fi)G(n)l  >  n 
where  in(fi)G(n)i  <  n  ' 


(20) 
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Depending  on  the  shape  of  the  sampling  pulse  go  (and  its  width  A)  and  on  the  noise  levels, 
the  frequency  bands  in  Ps  where  |G(f2)|  <  fi  may  be  large.  In  these  bands,  there  is  a  total 
loss  of  information  about  the  underlying  signal. 


Another  method  of  regularization  is  to  add  a  small  positive  constant  to  the  denom¬ 
inator  of  the  MNLS  solution  (15).  This  approach  is  popularly  known  as  the  Tikhonov 
regularization  method.  The  MNLS  solution  with  Tikhonov  regularization  is  given  by 

-1...  (21) 


xt  =  T*{TT*  +  fiB)-^yd 


In  the  frequency  domain  the  solution  is 


Ndi^r) 


li>0. 


(22) 


n(fi)|G(fi)P  +  /x  ^U{Q)\G{Q)\^  +  fi 

The  following  example  illustrates  the  two  types  of  regularization  schemes  and  their  effects 


on  the  resolution  limit- 


Example  3.1  Let  x{t)  be  bandlimited  to  [-IStt,  167r]  with  a  spectrum  as  shown  in  Figure  2. 
Let  the  sampling  function,  go{t)^  be  a  square-pulse  of  width  A  =  ^  and  let  r  =  ^ 

Figure  2  shows  a  plot  of  n(f2)|G(n)|.  With  Nd  equal  to  a  constant  p  =  10  and  a 
regularization  parameter  fi  =  25p,  the  two  regularized  solutions  Xi  and  Xr  are  illustrated 
in  Figure  2.  * 

In  the  absence  of  multiple  exposures  (with  r  =  A),  the  spectrum  of  X(fi)  could  have 
been  recovered  only  up  to  frequencies  in  [■^j  Because  of  multiple  interleaved  exposures 
(r  =  A/4),  we  are  able  to  extract  information  at  the  higher  frequencies  as  well.  From 
Figure  2  it  is  evident  that  Xr  is  zero  over  certain  frequency  bands.  It  is  also  clear  that 
when  the  regularization  parameter,  /i,  is  decreased,  the  size  of  these  missing  frequency 
bands  decreases;  however,  the  contribution  of  noise  to  the  solution  is  increased.  Thus 
the  parameter  fx  presents  a  trade-off  between  these  two  effects  and  must  be  chosen  as  a 
compromise  between  the  two.  The  zeroing  of  Xr  over  these  bands  (however  narrow)  may 
be  unacceptable  in  many  applications,  since  it  causes  oscillatory  behavior  (often  called 
ringing)  in  the  time  domain.  With  the  Tikhonov  regularization,  it  is  observed  that  the 
signal  component  of  the  MNLS  solution  is  approximate  in  the  entire  frequency  region. 
However,  there  are  no  missing  bands  and  the  time-domain  solution  is  smoother.  Here  as 
well,  choice  of  /x  is  critical  to  the  reconstruction  and  is  usually  made  based  on  the  SNR. 

Since  in  the  noise-corrupted  case,  the  signal  cannot  (in  general)  be  perfectly  recon¬ 
structed,  a  new  measure  of  resolution  is  necessary  that  allows  for  imperfect  reconstruction. 
We  develop  a  measure  of  resolution  based  on  the  maximum  tolerable  worst-case  error. 
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Figure  2:  Effects  of  the  spectral  truncation  and  Tikhonov  regularization  in  the  frequency 

domain.  - true  spectrum; - with  spectral  truncation;  —  •  —  with  Tikhonov  scheme; 

•  •  •  sampling  kernel  spectrum. 

It  is  obvious  that  the  magnitude  of  the  reconstruction  error  (even  after  normalization 
by  signal  norm)  depends  on  the  magnitude  of  the  measurement  noise  and  on  the  particular 
underlying  signal.  We  therefore  assume  that  the  SNR.  is  greater  than  a  given  positive 
constant  In  other  words,  we  assume  that  in  (4),  ||n(i|p  <  7/^||rx|p.  Therefore,  we  define 

resolution  as  follows. 

Definition  3.2  A  certain  reconstruction  algorithm  has  ^-resolution  of  5  or  better  if  the 
worst-case  normalized  reconstruction  error  (over  all  x  £  Bs  and  all  such  that  ||n(f||  < 
7/||rz||)  is  no  larger  than 

fcliE  <  £-  Vi  e  Bj;  Vn,  :  ||njf  <  ,=||ri||". 

Ikir 


Note  that  e^-resolution  is  defined  for  a  recovery  algorithm  and  not  for  the  recovery  problem. 

In  order  to  lower-bound  the  V-resolution  limit  for  a  particular  recovery  algorithm,  we 
require  tight  upper  bounds  for  the  worst-case  normalized  reconstruction  error.  As  we  shall 
see,  in  many  situations  the  upper  bound  might  trivially  be  1,  unless  the  feasible  set  of 
signals  is  restricted.  For  example,  let  5  =  :  |<j(n)|  <  /i}.  Consider  a  signal  x  &  Bs  such 
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that 


-(«)={; 

For  this  signal  it  is  clear  that  Xr(fi)  =  0  Vfi.  Thus,  if  a:  is  not  identically  zero,  the  worst-case 
normalized  error  in  this  case  is  1.  A  less  trivial  bound  on  the  worst-case  error  might  be  found 
if  (J  is  increased  or  if  we  restrict  the  set  of  feasible  signals  by  imposing  constraints  based  on 
known  signal  properties,  such  as  energy  bounds,  positivity,  bounds  on  the  derivatives,  or 
spectral  bounds.  We  will  attempt  to  obtain  this  bound  in  the  next  section  for  the  practical 
situation  in  which  only  finitely  many  measurements  are  available. 

4  Resolution  in  the  Practical  Situation 

In  this  section  we  analyze  resolution  limits  for  signal  recovery  problems  with  finitely  many 
noisy  observations.  Let  p  be  the  number  of  measurements  available.  The  measurement 
process  as  before  is  represented  by  the  linear  operator  T  :  Lj  (R) 

{Tx)i  =  (z,  flr,),  i  =  1, 2, . . . , p. 

The  Qi  could  be  uniformly  translated  versions  of  a  single  measurement  signal  pot  as  in  (7), 
but  it  is  not  necessary  for  this  analysis. 

The  reconstruction  problem  is  formulated  as  follows 

given  T,  S,  and  measurement  vector  j/j  €  C^,  find  x  E  Bg  such  that  Tx  =  yd- 

The  set  of  all  signals  in  Bg  satisfying  Tx  =  ydvsa  linear  variety,  V,  with  finite  codimension 
p.  Thus,  there  are  an  infinite  number  of  feasible  solutions  if  data-matching  is  the  only 
constraint.  The  min-max  optimal  solution,  which  is  also  the  minimum  norm  solution,  is 
given  by 

xmn  =  arflfmin||z|| 
x€V 

=  T*{TT*)-^yd. 

The  operator  TT*  is  now  simply  a  pxp  matrix  whose  ij*^  entry  is  {gj,  gi). 

Since  the  true  signal  can  be  any  member  of  V,  the  su premum  of  the  normalized  recon¬ 
struction  error,  >  is  1.  Thus  our  earlier  definition  of  e'- resolution  becomes  mean¬ 

ingless,  since  for  every  yd  and  passband  Pg  =  [^,  j],  we  can  find  a  signal  x  €  Bg  for  which 
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the  normalized  reconstruction  error  is  greater  than  the  tolerable  error.  Hence,  we  have  to 
restrict  the  set  of  admissible  signals  appropriately  to  bound  the  worst-case  error. 

However,  a  Unite  number  of  measurements  is  often  justified,  because  most  signals  en¬ 
countered  in  physical  systems  are  essentially  timelimited.  Accordingly,  we  restrict  our  at¬ 
tention  to  these  signals. 

Let  r  be  a  compact  set  (usually  a  union  of  possibly  discontiguous  closed  intervals)  in 
R.  Let  W  :  L-iiR)  ->  Li^R)  denote  the  windowing  operator  to  T, 


Wx{t)  = 


{ 


i(t)  t  €  r 

0  else 


(23) 


Since  {Wx,y)  =  {x,Wy)  'ix,y  ^  L2{R),Wvs  ^  self-adjoint  operator.  A  signal  is  said  to 
be  e-essentiaJly  time-limited  to  F  if  ||W^x|p  >  (1~^)II®IP'  ^«,<5(r')  denote  the  set  of 

signals  which  are  bandlimited  to  [■^j  j]  €-essentially  timelimited  to  F,  i.e.. 


G.,5(F)  =  {x  €  :  ||T^x||2  >  (1  -  e)||x||2}.  (24) 


The  set  ^(.^(F)  represents  most  signals  encountered  in  physical  imaging  and  information 
systems.  Hence,  we  state  the  following  definition  of  resolution. 


Definition  4.1  A  reconstruction  algorithm  on  concentration  window  F  and  concentration 
factor  1-6  with  SNR  >  ^  will  be  said  to  have  c'-resolution  of  S  or  better  if  Vx  6  G«,5(F) 


-  and  Vnj  s.t.  ||n<i||  <  »7lir®ll> 


lla:  -  x|l^ 

IkIP 


<  e'. 


The  set  Ge,5(F)  has  several  interesting  properties  which  can  be  exploited  to  determine  res¬ 
olution  limits.  Many  of  these  properties  are  characterized  by  an  orthonormal  family  of 
functions  called  the  Prolate  Spheroidal  Wave  Functions  (PSWFs), 

sociated  eigenvalues,  {A,}gi  [16],  [17].  The  PSWFs  corresponding  to  F  and  Pj  are  solutions 
to  the  following  integral  equation. 

Properties  of  the  PSWF  and  the  associated  eigenvalues  have  been  widely  studied.  The 
following  three  properties  of  the  PSWF  are  essential  to  the  treatment  of  resolution  limits 
with  finite  data: 
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1.  <j>i’s  are  bi-orthogonaJ  functions,  i.e.,  for  i  ^  j 


{W(f>i,W<i>j)  =  0 


2.  Ai  >  A2  >  As  •  •  •  >  0, 

3.  form  an  orthonormal  basis  for  Bs^ 

With  Ge,j(r)  as.  the  feasible  set,  the  recovery  problem  becomes 

given  T,  5  and  the  measurements  yj  G  C’’,  find  x  €  Ge,5(r)  such  that  Tx  =  yd- 

Since  the  resolution  limit  of  a  recovery  algorithm  is  based  on  the  worst-case  relative  error, 
our  objective  is  to  find  an  algorithm  that  will  minimize  the  worst-case  relative  error. 

4.1  Lower  Dimensional  Approximation 

Consider  any  signal  x  €  Ge,5(r)  C  Bs-  Since  the  form  an  orthonormal  basis  for  Bs 

from  property  (3),  we  can  express  x  as 

00 

I  =  ^ 
t=i 

To  recover  x  6  Ge,5(r),  in  general,  we  have  to  determine  an  infinite  number  of  a,-,  from  a 
finite  number  (p)  of  observations.  However,  every  x  €  C?e,5(r)  satisfies  >  1  “  a-iid 

hence  we  have  the  following  condition  on  the  coefficients  a,-, 

^o-Ai  >  1  -  €. 

•=i 

This  condition  suggests  that  we  may  be  able  to  restrict  our  reconstructions  to  a  finite¬ 
dimensional  subspace  of  Bs  (of  dimension  less  than  p)  and  still  obtain  low  error  recon¬ 
structions.  The  next  theorem,  a  well-known  result,  shows  that  G£,i(r)  is  essentially  finite 
dimensional  and  that  the  PSWF  optimally  approximate  Ge,j(r)  [16],  [17]. 

Theorem  4.1  (Landau-Pollak-Slepian)  For  any  positive  integer  N,  amongst  aJl  N- 

dimensional  spaces  Sn,  the  space  spanned  by  the  first  N  PSWF,  Sf^  =  span{<l)i  (Pn}, 

is  optimal  for  Ge,j(r)  in  that  it  minimizes 

max  min  — jj — jj — .  (zoj 

*€G«,i{r}  *€Sjv  IfII 
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Moreover,  the  worst-case  relative  error  can  be  expressed  in  terms  of  the  eigenvalues  associ¬ 
ated  with  the  PSWF  as 


1  f  1  ^N+l  >  1  —  €  >  0 


(27) 


E{Sff)  is  also  called  the  N -width  of  the  set  G«,5(r)  or  the  intrinsic  error  at  dimension  N . 


Thus,  for  a  fixed  dimension,  r,  the  subspace  spanned  by  the  minimizes  the  worst- 

case  relative  error.  Moreover  the  worst-case  relative  error  E{S^)  decreases  with  r.  Thus  it 
would  seem  that,  given  p  observations,  the  best  choice  of  a  lower  dimensional  subspace  to 
approximate  (Je,5(r)  would  be  =  span{^i,  •  •  • ,  <i>p}~  This  choice  would  lead  to  p  equations 
in  p  unknowns.  Unfortunately,  thep  parameters  required  to  describe  the  reconstruction  from 
Sf  cannot  be  determined  exactly  from  the  observations  yd  for  two  reasons.  First,  yd  are 
noise-corrupted  in  practice.  Second,  the  observations  yd  are  linearly  related  to  i  € 
and  not  to  the  projection  Xp  of  x  onto  5^.  Thus,  an  additional  error  will  be  incurred  in 
determining  the  parameters  that  describe  the  lower  dimensional  approximate.  We  next 
derive  an  expression  for  this  error  and  its  worst-case  value  H(r)  for  a  fixed  dimension  r .  We 
suggest  choosing  r  to  minimize  E{Sf)  +  H(r). 

4.2  Worst-Case  Error  Analysis  for  Subspace  Selection 

Consider  the  reconstruction  based  on  an  r-dimensional  approximation  of  G«,5(r),  where 
r  <p,  and  let  Xr  be  the  projection  of  x  onto  Sf.  Then, 

Xr  =  Y^Oli<f>i  (28) 

«=i 

OO 

and  the  approximation  error  ®  ~  ®r  =  lii  section,  we  study  the  effect 

i=r+l 

of  measurement  noise  and  approximation  error  Cr  on  the  estimate  of  from  measure¬ 

ments  yd-  The  measurements  yd  are  linearly  related  to  x  and  corrupted  by  noise  Ud- 

yd  =  Tx  +  nd 

=  Txr -i-Tcr -ir  nd 
r 

=  ^2  +  (TCr  -h  nd) 

i=l 

=  Ara’’  -\-{Ter-\-nd),  (29) 
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where  Ar  is  a  p  x  r  matrix  with  Ar{i,j)  =  {<l>j,gi),  We  will  assume 

that  the  columns  of  Ar  are  linearly  independent.  If  they  are  not,  the  worst-case  error  is 
unbounded.  The  LS  solution  of  Ti  =  y,i,  i  €  Sf,  is  determined  from  the  MNLS  solution 
of  Aro:'’  =  yd,  which  is 

d’’  =  A^yd 

=  o’"  +  aI  {Ter +  nd), 

where  aJ  =  (A^Ar)“^A^.  Thus  the  reconstruction  is  given  by 

Xr  =  ^Qi<f>i  (30) 

j=i 

r 

and  Cr  =  a:r  -  ir  =  is  the  additional  error  incurred  in  determining  the  a* 

i=l 

parameters.  The  error  term,  Cn  has  contributions  from  both  Ter  ‘^nd  measurement  noise 
rid- 

Thus  the  total  reconstruction  error  is 

=  \\X  -  Xr  +  Xr  -  Xrf 

=  Ikr+Crf 

=  ||er||^  +  IlCriP  since  Cr  X  sf  and  obviously  Cr  €  Sf. 

A  pictorial  representation  of  these  error  terms  is  given  in  Figure  3.  We  seek  an  upper  bound 


Figure  3:  Error  terms. 
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on  IK.  From  Theorem  4.1,  it  is  clear  that  V®  6  (j«,i(r), 


^N+l  >  1  >  0 

Aiv+i  <  1  -  €  <  Ai 


All  that  remains  to  be  done  is  to  upper-bound  |^j|!r »  This  error  term  can  be  interpreted 
as  a  penalty  on  complexity.  The  complexity  in  this  case  is  the  dimension  of  the  subspace. 

Recall  that  Cr  =  -  0£i)<f>i  and  that  d*"  -  =  i4;(rer  +  n^).  Hence, 

tsl 


IICIK  =  IK-alK 

=  WAtiTer  +  nM^ 

<  UlTerf  +  WAUdf  +  2|lAtre,||  ||4«<ili,  (31) 

with  equality  achieved  when  aJtij  is  collinear  with  A^Ter-  If  a  lower  bound  (;^)  on  SNR 
is  known,  then 

<  ^WTxf 

<  >)V„„(rr)  iixip, 

with  equality  achieved  when  x  is  the  singular  vector  of  T  that  corresponds  to  its  maximum 
singular  value  y/ CTmax  (TT*) .  Now  all  that  remains  to  be  done  is  to  obtain  a  tight  upper 

bound  on  over  x  €  Gc^^(r).  A  loose  upper  bound  may  be  obtained  by  using 

II^Terll  <  but  we  can  achieve  a  better  bound  because  Cr  has  less  than 

of  its  energy  in  F.  Mathematically,  our  objective  is  to  find 

sup  II  f;  aiAtT<f>i\\\  (32) 

This  is  a  nonlinear  infinite  programming  problem,  and  in  general  we  can  seek  only  approx¬ 
imate  solutions.  We  devote  the  remaining  of  this  section  to  that  task.  We  first  show  that 
under  mild  conditions,  the  problem  can  be  approximated  by  a  finite  variable  nonlinear  pro¬ 
gramming  problem.  We  then  apply  well-known  techniques  to  solve  the  problem.  We  start 
with  the  following  definition. 
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Definition  4.2  An  observation  made  with  a  measurement  functional,  g,  is  in  the  concen¬ 
tration  window,  r,  if  support(g)  C  F 


Lemma  4.1  [18]  The  eigenvalues  {A,}  associated  with  the  PSWF  corresponding  to  a  pass- 
band  [^,  j]  and  a  compact  concentration  window  F  are  absolutely  summable  and 


(33) 


Lemma  4.2  Let  x  €  Ge,5(F)  with  x  =  °‘k<f>k  and  ||x|p  =  al  =  1.  Let  ejv  =  ^  ak(l>k 

*=1  k=N+l 

If  all  the  measurements  are  taken  in  the  concentration  window,  F,  i.e.,  every  gi  has  a  sup¬ 
port  inside  the  concentration  region  F,  then  there  exists  a  finite  number  M  such  that  for 
every  positive  integer  N, 

||Tejv|P<M(  f;  Afc),  - 
k=N+l 

and  lim  \\TeN\f  =  0.  In  fact,  ||TT*1|  will  serve  as  M.  ■ 

N-¥oo 

Thus  as  N  becomes  large  the  contribution  of  Tejv  becomes  negligible  in  the  error. 


Theorem  4.2  Let  7  >  0  and  the  rest  of  the  assumptions  be  as  stated  in  Lemma  4.2.  There 
exists  a  positive  integer  N  independent  of  the  choice  of  x  E  C?«,5(F)  such  that 


Furthermore, 


N 


Il4re,||  -  II  S  “i4r*||  <  T- 

i=r+l 


(34) 


sup 

x€Ge,j(F) 


UkerW 

Ikll 


=  7  + 


SUP  II  X]  “i^r^’^tll-  (35) 


Proofs  of  Lemma  4.2  and  Theorem  4.2  can  be  found  in  the  appendix. 

Theorem  4.2  shows  that  the  solution  to  the  infinite  programming  problem  in  (32)  can 
be  approximated  closely  by  a  suflSciently  large  finite  variable  problem.  Let  N ,  sufficiently 

large,  be  chosen,  let  Hjsf  =  ^2?  •  •  be  a  p  ><  iV  matrix  with  columns  hi  given  by 
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^  ^  r  0  l<»<r 

\  ^  1  ^  i  ^ 


and  let  An  —  diag([Ai,A2, •••.Aat]).  Then  a  bound  on  \\A^Ter\\  can  be  obtained  by  com¬ 
puting 

nv=  max  ||Ifiva^||  (37) 

||A*aJV||2-l-« 

A  popular  method  for  solving  this  nonlinearly  constrained  quadratic  program  is  the  sequen¬ 
tial  quadratic  programming  method,  [13],  [19]. 

Let  6,,  =  (7 -I- and  define  the  noise  contribution  factor  .  Then 

^m«n 

IMrre.(p  <  6r|l®l|*  and  from  (31), 

IKrII’  <  Il4re.||“  +  ||4»dll“  +  2|l4r€.||||4n4| 

(h  \  )_2  ,  o  f^masiTT*)  ^  ^  2 

=  ibr  +  pW  +  2^/b^)\\xf 

=  =('•)ll«ll^  f38) 


where  S(r)  _  br+Prif +2\/brP^r]^.  Hence,  we  have  the  following  procedure  for  determining 
S(r),  using  the  PSWF,  <bi,  the  corresponding  eigenvalues  A,-,  and  an  acceptable  value  of  7 

For  each  candidate  dimension  r, 

•  construct  the  p  x  r  matrix  Ar  =  [T(j>iT<f>2'  •'T<l>r]  and  compute  its  smallest  singular 
value,  (Tmin. 

•  Find  the  smallest  integer  N  for  which 


S  Jr  ft  TT* 


•  Construct  Hn  as  in  (36),  and  let  An  =  diag([Ai,A2,-**,AAr]).  Solve  (37)  by  the 
sequential  quadratic  programming  method  and  let  br  =  (7-f  mr)^. 

=(r)  =  ir  +  Pr’7^  +  2y^!p^=  (v^-|-pr77)^  (39) 
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Thus  a  bound  on  the  worst-case  normalized  error,  ©(r),  can  be  obtained  by  the  sum  of 
the  intrinsic  error  and  a  bound  on  the  worst-case  ||Cr|P) 

0(r)  =  min(l,  E{Sl)  +  H(r)).  (40) 

Thus,  to  ascertain  whether  a  resolution  of  5  can  be  achieved  given  a  set  of  p  measure¬ 
ments,  we  first  compute  the  bound  on  the  worst-case  normalized  error  bounds  for  each 
dimension  r  ranging  from  1  through  p  using  the  PSWF  corresponding  to  Pj  =  [ir,  f  ]  and 
we  determine  the  dimension  r*  which  gives  the  smallest  error.  If  the  worst-case  error  for  this 
dimension  is  below  the  allowable  error,  we  can  claim  that  a  resolution  of  5  can  be  achieved 
with  the  given  set  of  measurements  and  noise  level. 

Remarks: 

1.  As  a  consequence  of  the  above  analysis  we  have  a  new  algorithm  for  signal  recovery 
based  on  dimension  reduction  guided  by  the  bound  on  the  worst-case  reconstruction 
error. 

2.  The  worst-case  error-bound  given  by  (40)  does  not  depend  on  the  data  yd-  It  depends 
only  on  the  sampling  functions,  gi,  the  bandwidth,  the  sample  spacing,  r,  the 
noise  level,  g,  and  the  choice  of  the  dimension,  r.  Thus,  the  selection  of  the  dimension 
and  the  determination  of  resolution  can  be  made  (off-line)  before  the  measurements 
are  taken. 

3.  Our  analysis  and  definition  of  resolution  are  based  on  worst-case  errors  in  a  determin¬ 
istic  framework.  Therefore,  in  general  the  reconstruction  error  can  be  expected  to  be 
lower  than  the  predicted  value. 

4.  The  analysis  holds  true  for  all  sampling  patterns.  Hence  gi,  the  measurement 
function,  does  not  have  to  be  a  shifted  version  of  a  single  measurement  function  go- 
The  only  restriction  is  that  the  support  of  each  g,  lies  inside  the  concentration  window, 

r. 

5.  We  have  assumed  essential  timelimitedness  and  strict  bandlimitedness  in  our  treat¬ 
ment,  which  can  be  easily  changed  to  essential  timelimitedness  to  F  and  essential 
bandlimitedness  to  P.  The  PSWF  will  still  be  the  optimal  sequences,  [20]  and  all  the 
results  will  still  hold  true,  with  minor  modifications. 

■  6.  The  PSWF  have  been  studied  in  the  classical  setting  of  ID  signals  with  lowpass 
passband  and  contiguous  concentration  intervals.  The  three  relevant  properties  of  the 
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PSWF  and  the  dimensionality  theorem  can  be  generalized  to  the  more  general  setting 
of  mD  signals  with  P  and  T  discontiguous.  They  follow  by  expressing  the  integral 
equation,  (25)  as  a  linear  operator  equation,  as  follows  [21],  [22] 

WBWi^i  =  Xi-ipi  ui) 

Xi(f>i  =  BWi)i  '  ^  ’ 

Recognizing  that  the  are  the  eigenvectors  of  a  positive  semidefinite  operator  with 

eigenvalues  A,-,  the  following  three  properties  of  the  PSWF  which  are  essential  to  this 

analysis  become  obvious: 

•  the  <f>i's  are  biorthogonal  functions,  i.e.,  for  i  ^  j 

{W<f>i,W(l>j)  =  0 


•  Ai  >  Aa  >  As  •  •  •  >  0. 

•  forms  an  orthonormal  basis  for  Bp,  the  space  of  finite-energy  signals 
bandlimited  to  P. 

Since  the  dimensionality  theorem  uses  only  these  three  properties,  it  can  be  generalized 
as  well. 


4.3  Implementation  Details 

In  practice,  all  integrals  have  to  be  computed  using  numerical  approximation  methods.  We 
use  simple  summation  after  discretizing  functions  on  a  fine  grid.  With  discretization,  the 
PSWF  computation  of  (41)  becomes  an  eigenvalue-eigenvector  computation  followed  by 
low-pass  filtering  [23].  Constructing  the  entries  of  the  matrices  Ar,  Hn  and  TT*  requires 
computation  of  integrals  and  is  also  achieved  by  discretization  and  summation. 

In  the  2D  case,  if  the  passband  is  square,  the  kernel  in  the  integral  equation,  (25), 
describing  the  PSWF  becomes  separable.  If  in  addition  the  concentration  region  is  also 
square,  the  double  integral  in  (25)  becomes 


IJik^k{t,3)=  /  /  - ^ ^ - p - ;-^^k{n,T2)dTidT2,  (42) 

^  J-T/2J-T/2  ff(i-n)  Trit-T2) 


where  is  the  2D  PSWF  with  associated  eigenvalue  fik .  If  the  <f>i  are  the  ID  PSWF  cor¬ 
responding  to  5  and  T  =  [-r/2,T/2],  then  it  is  clear  that  for  each  i,j,  $(i,s)  =  (^,(t)^j(s) 
satisfies  (42)  with  (x  =  A^Aj.  Thus  after  discretization,  if  we  represent  each  2D  PSWF 
as  a  long  ID  vector  of  columns  stacked  one  below  the  other,  then  the  2D  PSWFs  can  be 
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computed  as  Kronecker  products  of  corresponding  ID  PSWFs  taken  two  at  a  time  with 
eigenvalues  equal  to  the  product  of  the  corresponding  eigenvalues. 

We  next  illustrate  these  results  and  the  computational  details  in  two  examples.  The 
first  example  is  in  a  ID  setting.  The  second  example  is  in  a  2D  setting  to  illustrate  the 
applicability  of  the  results  to  multiple  dimensions. 

Example  4.1  (ID  setting)  Consider  reconstruction  of  ID  signals  that  are  bandlimited 
to  [-47r,47r]  and  have  at  least  99.5%  of  their  energy  concentrated  in  [-2.0, 2.0],  from  19 
measurements  taken  with  shifted  unit  rectangular  pulses  of  width  ^  and  interpulse  distance 
(stagger)  Thus  8,  T,  and  €  are  i,  [-2.0, 2.0]  and  0.005  respectively  in  the  definition  of 
Ge,5(r),  while  p  and  r  are  19  and  The  sampling  functions  are 

!7Jfe(0=5o(*-^)  fc  =  l,2,...19. 


where 


9o{t)  = 


{ 


1 

0 


0  <  f  <  0.5 
else 


Note  that  the  width  of  the  sampling  pulses  is  i  =  28:  Thus  the  Rayleigh  resolution  limit 
is  28.  Let  the  error  tolerance  be  10%.  We  will  see  that,  at  an  SNR  of  40dB,  a  resolution 
limit  of  r  =  i  (better  than  the  Rayleigh  limit)  can  be  achieved  by  the  proposed  algorithm. 

The  nonlinear  programming  problem  of  (37)  is  solved  by  using  the  sequential  quadratic 
programming  method.  The  intrinsic  error  E{Sf),  the  parameter  estimation  error  without 
noise  br,  the  noise  contribution  factor  pr,  are  computed  for  each  dimension,  r,  for  which 
E{Sf)  <  1.  The  values  are  tabulated  in  Table  4.1  along  with  0(r)  =  E{Sf)  +  Pr'^Y 

for  SNR  of  32  dB  and  40  dB.  From  the  table  we  observe  that,  with  a  40dB  SNR,  the 
optimal  dimension  for  this  signal  recovery  problem  is  15,  and  the  worst-case  normalized 
error  is  bounded  above  by  0.0531.  In  fact,  using  (39,  40),  we  can  show  that  with  a  10% 
error  allowance,  a  resolution  of  at  least  r  —  0.25  can  be  achieved  by  the  above  algorithm 
whenever  the  SNR  is  greater  than  32dB.  As  remarked  earlier,  all  these  computations  can 
be  done  off-line,  since  they  do  not  depend  on  the  actual  observed  data. 

We  now  test  the  reconstruction  algorithm  using  the  precomputed  optimal  dimension  of 
15  on  a  specific  signal. 


x{t)  =  ^ 


_  ^sin(0.47rt)  Y 
Q.Airt  ) 


-1-0.2 


sin(O.lfft) 

O.lTTt 


cos(3.57rt). 
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Table  1:  Intrinsic  error,  noise  factor  and  total  error  versus  dimension. 


■ 

~wsfr 

br 

noise-factor 

Pr 

e(r) 

SNR=40dB 

e(r) 

SNR=32dB 

0.2203 

0.1905 

4.3953 

0.4511 

0.5188 

HI 

0.0379 

0.0430 

5.4552 

0.1065 

0.1561 

15 

0.0100 

0.0189 

7.0128 

0.0531 

0.1078 

16 

0.0046 

0.0338 

9.4899 

0.0823 

0.1819 

17 

0.0034 

0.0788 

,13.2808 

0.1744 

0.3788 

18 

0.0032 

0.1103 

17.6118 

0.2615 

0.5998 

19 

0.0031 

0.1124 

27.8416 

0.3797 

1.0000 

A  plot  of  the  signal  is  shown  in  Figure  4  (a).  The  highest  frequency  in  x{t)  is  3.6  radians. 
Note  that  a?(t)  is  a  low  frequency  signal  with  a  low  energy,  high-frequency  ripple.  The 
frequencies  are  selected  in  order  that  the  low-frequency  component  falls  below  the  Rayleigh 
limit  and,  hence,  is  captured  by  the  observations,  whereas  the  high-frequency  ripple  is  much 
above  the  Rayleigh  limit  and  thus  is  not  seen  in  the  observations  (Figure  4  (b)).  A  high 
resolution  reconstruction  should  resolve  the  high  frequency  ripple,  also.  Since  x{t)  has  99.69 
%  of  its  energy  inside  the  concentration  interval  [—2.0, 2.0],  i.e.,  e  =  0.0031,  it  belongs  to 
the  set  G4,5(r)  considered  in  this  example.  We  take  19  observations  in  the  concentration 
window  with  shifted  versions  of  the  sampling  function  described  by  (4.1).  The  observations 
are  shown  in  Figure  4(b).  Note  that  the  high-frequency  ripple  is  completely  lost  in  the 
observations. 

The  reconstruction,  x,  is  computed  using  the  algorithm  with  precomputed  dimension 
of  15.  It  is  also  depicted  in  Figure  4  (a).  The  normalized  error  for  this  reconstruction, 

(Hj|~|ll)2,  is  computed  to  be  0.0068,  which  is  much  less  than  the  worst-case  error  bound  of 
0.0531.  ■ 

We  now  consider  an  example  in  two  dimensions. 

Example  4.2  (2D  setting)  Consider  the  reconstruction  problem  for  2D  signals  that  are 
bandlimited  to  a  square  passband  Ps  —  [— 27r,2jr]  x  [— 27r,  27r]  and  have  at  least  99%  of 
their  energy  inside  the  square  region  [—2.5, 2.5]  x  [—2.5, 2.5]  from  121  measurements.  The 
measurements  were  taken  using  pulse  sampling  functions  with  a  square  region  of  support  of 
width  1  and  an  intersampling  distance  of  5  in  each  direction.  Since  the  sample  spacing  (5)  is 
smaller  than  the  sampling  kernels’  width  (=1),  this  sampling  strategy  requires  interleaving 
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Figure  4:  Reconstruction  of  the  test  signal  using  the  optimal  dimension  determined  by  the 

worst-case  error  analysis.  The  SNR  is  40db.  (a) - signal,  x(t); . reconstruction,  £ (t) 

(b)  +  Observations. 

of  multiple  (four)  staggered  exposures  as  in  Figure  1.  The  sampling  functions  are 


/  k 

9l,ki^,s)  =  go{t- -,s- -)  /,*  =  l,2,...ll. 


where 


5o(i,s) 


1  0<«<l,0<s<l 
0  else 


With  this  sampling  kernel,  the  Rayleigh  limit  will  be  1  in  each  dimension.  Let  measurements 
be  taken  from  the  concentration  window  with  t}  =  This  value  of  tj  corresponds  to  an 
SNR  of  40  dB. 
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The  nonlinear  programming  problem  of  (37)  is  solved  using  the  sequential  quadratic 
programming  method,  and  E{Sf),  H(r),  and  0(r)  are  computed  for  each  dimension  r  for 
which  E{Sf)  is  less  than  1.  These  values  are  plotted  in  Figure  5.  From  the  figure  we 


Figure  5:  Intrinsic  error,  E{r) - ;  parameter  approximation  error,  E(r) - ;  total  error. 


observe  that  the  optimal  dimension  for  this  signal  recovery  problem  with  a  40dB  SNR  is 
any  dimension  between  81  and  89.  The  worst-case  normalized  error  in  the  noisy  situation 
(40dB  SNR)  is  0.1326.  Thus,  if  an  error  of  13.5%  is  tolerable,  then  a  resolution  of  r  =  0.25 
can  be  achieved  by  the  above  algorithm. 

Now  consider  a  specific  signal  in  Ge,i(r):  x{t,s)  =  h{t)h{s)  with 


h{t)  =  0.25 


/  stn(0.2fft)'\^ 
V  0.27ri  J 


+  0.05 


sin{0.2Trt) 

0.2irt 


cos(1.757rf). 


A  meshplot  of  the  x{t,  s)  is  shown  in  Figure  7.  As  in  the  ID  case,  the  signal  is  comprised  of  a 
low-frequency  component  with  a  high-frequency  ripple.  The  observations  are  four  exposures 
with  stagger,  and  they  are  as  shown  in  Figure  6.  The  interleaved  observations  are  also  shown 
in  Figure  6.  The  reconstruction  using  dimension  81  is  shown  in  Figure  7.  The  relative  error 
is  computed  to  be  0.0179,  which  is  much  below  the  computed  upper  bound.  ■ 
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5  Conclusions 


A  resolution  analysis  for  signal  recovery  from  finitely  many  discrete,  noise-corrupted,  linear 
measurements  is  presented.  A  new  measure  for  resolution  is  introduced,  which  is  more 
appropriate  than  the  Rayleigh  resolution  limit  for  current  signal  recovery  algorithms.  This 
resolution  measure  is  based  on  a  prescribed  tolerance  of  relative  error  in  the  reconstruction, 
and  unlike  previous  definitions  is  able  to  bring  out  the  extent  to  which  time  or  spatial  domain 
features  can  be  recovered  by  an  algorithm.  The  computation  of  resolution  limits  reduces 
to  the  computation  of  the  worst-case  relative  error  in  the  recovered  signal.  By  suitably 
constraining  the  class  of  feasible  signals,  the  worst-case  error  is  expressed  as  the  solution  of 
a  finite-variable  nonlinear  program.  The  analysis  and  examples  show  that  details  finer  than 
the  Rayleigh  resolution  limit  can  be  recovered  by  simple  linear  processing  even  in  practical 
situations  with  finite,  noise-corrupted  data.  In  the  process,  we  derive  an  algorithm  for  high 
resolution  reconstruction  (from  linear  observations)  and  show  how  one  can  precompute 
worst-case  error  bounds  and  the  resolution  limit  for  the  algorithm. 

Appendix 

Proof  of  Lemma  4.1:  Since  all  the  measurements  are  taken  inside  the  concentration 
region,  Tx  =  TWx,  Vx  €  C?e,5(r).  Therefore, 

llrewf  =  II  f;  c.tT4,tf 

=  \\T  f;  akW<i>kf 

<  lirr-||  f;  al\\w<t>k\\^ 

k=N+l 

<  M  ^  Xk  where  M  =  UrT’^H 

k=N+l 

Lemma  4.1  and  property  (2)  of  the  PSWF  together  show  that  =  0*  • 

Proof  of  Theorem  4.2:  Let  (Tjnin  ^  0  he  the  smallest  singular  value  of  Ar>  By  Lemma  4.1 
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it  is  possible  to  choose  N  such  that  ^  Xk  <  Then, 

k=N^l  ^  " 


N  oo 

\\4Ter\\  =  II  E  <^i4T<l>i+  E  «*-4^<^ill 

i=N’+l 


<  II  x;  «i4r0ill  +  ll42’eiv|l 

i=r+l 


<  II  f;  aiAt^’<^ill  +  ll4lll|3’eiv|| 

»=r+l 

<  II  Y'  aiAjr<^;||  +  -^.  M  ^  (by  Lemma  4.2) 

t=r+l  ^  i=N+l 

<  II  x;  «.4r<Ai||  +  7 

J — »_1_1 


Next, 


sup  =  sup  II  x:  <^4T4>i\ 


<  sup 

<  sup 
X)ili 


II  x;  «fAtr?i.||  +  T  (from  43) 

i=:r+l 

II  £  «i4r«H-7 


^  + 

=  sup  II  Y  «i4r0.il  +  7  (■ 

Er.i“?=i  ‘=’-+' 

Inequality  (44)  follows  from  the  fact  that  the  set  5i  =  {a  :  <*t  =  1;  ESi  (1  — 

e}  is  contained  in  the  set  52  =  {a  :  Ei^i  —  liEi^i  ~  ‘^«)  —  hence 

supremum  over  S2  is  greater  than  the  su premum  over  Si.  Equality  (45)  follows  beca 

II  Eilr+i  «t4r^t||  is  a  convex  function  in  a,  and  a  convex  function  maximized  ove 
convex  set  achieves  its  maximum  at  the  boundary  of  the  convex  set. 
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Figure  7;  Signal,  x{t,s),  and  reconstruction,  x(t,s). 
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(e) 


Figure  6:  Observations:  (a),  (b),  (c),  (d)  Four 


observations. 


